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This  Symposium  has  dealt  with  provable  algorithms  for  finding  zeros 
of  general  polynomials,  with  the  tacit  assumption  that  the  processes 
would  be  implemented  on  an  ideal  computer  system  capabl®  of  exact  arith¬ 
metic  operations.  In  contrast,  I  should  like  to  point  out  the  near 
absence  of  algorithms  to  solve  even  a  quadratic  equaticn  in  a  satisfactory 
way  on  actually  used  digital  computer  systems.  The  difficulties  ai-e 
partly  caused  by  round-off  error  in  floating-point  arithmetic,  but  much 
more  by  the  ever-present  possibility  of  overflow  or  underflow  (defined 
below) . 

This  note  presents  specifications  for  a  satisfactory  quadratic 
equation  solver  suggested  by  Professor  W.  Kahan  of  the  University  of 
Toronto  in  lectures  at  Stanford  University  in  1966.  The  general  level 
of  performance  is  Kahan's,  but  the  details  are  mine. 

Consider  the  following  set  F  =  F(3,  t,  m,  M)  of  normalized  floating' 
point  numbers.  This  uses  a  niumber  base  3  (bases  2,  8,  10,  and  I6 
are  in  use)  and  a  prescribed  number  t  of  significant  digits.  There 
are  two  limiting  integer  exponents  m  and  M  .  The  set  F  contains 
precisely 


^Tb.is  is  a  brief  written  contribution  to  Symposium  on  Constructive 
Aspects  of  the  Fundamental  Theorem  of  Algebra,  held  at  the  IBH  Research 
Lcbcratories,  .Rueschlikon,  Switzerland,  5-7  June  196?.  The  preparation 
or  this  note  and  the  author's  attendance  at  the  Symposium  were  sponsored 
i  V  the  !J.  s.  Of  fice  of  Naval  Research  under  project  NR  044  211. 
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numbers.  One  of  these  is  0  .  Each  other  number  in  F  has  a  unique 
representation 

(2)  +  N  X  3®, 

where  the  sign  is  +  or  where  the  integer  exponent  e  aatiafies 
the  inequality 

(3)  m  <  e  <  M, 

and  where  the  integer  significand  N  satisfies  the  normalization  condition 

(4)  0^"^  <  N  <  0^  -  1  . 

Frequently  0  is  given  the  unique  computer  representation 

(5)  +0X0®. 

See  [l]  for  a  discussion  of  this  floating-point  number  system  in  a 

slightly  different  notation. 

* 

We  choose  F  to  be  a  certain  subset  of  F  consisting  of  numbers 
not  too  close  to  overflow  or  underflow. 

Definition.  To  be  definite  (and  somewhat  arbitrary),  let 
*  * 

F  =  F  (0,  t,  m,  M)  be  the  set  consisting  of  C  and  all  numbers  (2) 
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subject  to  (4)  and  also  to 


(3*)  m  +  l<e<M-l  . 

Definitions.  A  real  number  x  is  said  to  be  in  the  range  of  F 
if  either  x  =  0  or 

(6)  <  |x|  <  (a^  -  . 

* 

A  complex  number  z  is  iji  the  range  of  F  if  both  Re(2)  and 
Im(z)  are  real  numbers  in  the  range  of  F  . 

One  similarly  defines  the  expression  the  range  of  F  for  real 
and  complex  numbers. 

The  statement  (in  computer  Jargon)  that  a  real  number  y  suffers 
from  over-  or  underflow  is  equivalent  to  saying  that  y  is  not  in  the 
range  of  F  . 

In  terms  of  these  concepts  I  now  give  specifications  in  the  form 
of  a  commented  heading  in  Algol  60  of  what  I  consider  to  be  a  satisfactory 
quadratic-equation  solver  for  use  with  a  processor  of  floating-point  num¬ 
bers  in  F(t,  t,  m,  M)  . 

procedure  QUADPvATIC  (a,  b,  c,  xl,  yl,  x2,  y2,  error) ; 

value  a,  b,  c;  real  a,  b,  c,  xl,  yl,  x2,  y2;  switch  error; 

2 

comment  We  are  solving  the  equation  az  +  bz  +  c  e  o,  for  arbitrary 

* 

input  parameters  a,  b,  c  in  F  .  Where  values  of  the  output 
parameters  are  not  specified,  they  are  irrelevant. 
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If  a  =  b  =  c  =  0,  exit  to  error[l],  since  all  complex  num¬ 
bers  z  satisfy  the  equation. 

If  a  =  b  =  0  and  c  /  0,  exit  to  error[2],  since  no  z 
satisfies  the  equation. 

Otherwise,  let  and  be  the  exact  roots  of  the  equation, 

numbered  so  that  [z^^l  <  \z^\  .  (If  a  =  0,  let  ^2  “  * 

Whenever  is  in  the  range  of  F  ,  set  xl  and  yl  to 

nujribers  in  F  close  (defined  below)  to  the  real  resp.  imaginary 
part  of  . 

Whenever  z^  is  in  the  range  of  F  ,  set  x2  and  y2  to 

numbers  in  F  close  to  the  real  resp.  imaginary  part  of  Zg  . 

Let  =  xl  +  i  X  (yl)  .  We  require  that  ■  0  (if 
Zj^  *  0  ),  and  otherwise  (again  being  somewhat  arbitrary)  that 
-  z^l  <  p  +  1  units  in  the  least-significant  digit  of  the 
significand  of  the  floating-point  representation  (2)  of  max  (|xl|, 
|yl|)  .  To  repeat  this  requirement  in  symbols,  if  entierfx] 
denotes  the  greatest  integer  <  x,  we  demand  that 

-  z^l  <  (P  +  1)  X 

where  e  =  entier[logamax  (|xl|,  |yl|)  -  t]  . 

p 

We  make  a  corresponding  requirement  of  ^2  =  x2  +  i  X  (y2)  . 

If  Zj^  is  not  in  the  range  of  F,  but  Zg  is  in  the  range 
of  F  ,  set  x2,  y2  as  above  and  exit  to  error[5]  . 

If  Zg  is  not  in  the  range  of  F  (including  the  case 

\  * 

Zg  =•  ;,  but  Zy  is  in  the  range  of  F  ,  set  xi,  yl  as  above 

and  exit  to  error[4]  . 
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If  neither  nor  is  in  the  range  of  F,  exit  to 
errorf5]  • 

If  a  root  is  in  the  range  of  F  but  not  in  the  range  of 

* 

F  ,  we  permit  either  an  indication  of  over-  or  underflow  via  the 
appropriate  exit  to  error,  or  a  determination  of  z^  with  an 
accuracy  satisfying  (?)  above. 

The  procedure  QUADRATIC  should  make  no  unnecessary  use  of 
multiple-precision  computation,  but  computation  with  2t  significant 
digits  is  essential  at  one  part  of  the  procedure,  to  achieve  the 
accuracy  (?).  End  of  comment; 

The  main  source  of  practical  difficulty  in  writing  the  procedure 
QUADRATIC  is  the  possibility  of  over-  or  underflow  in  many  places.  What 
is  actually  programmed  depends  crucially  on  what  the  computing  system 
does  in  case  of  over-  or  underflow.  An  ideal  system  permits  the  user's 
program  to  regain  control  of  the  algorithm,  if  the  user  wishes,  with  the 
ability  to  interrogate  Boolean  variables  to  learn  whether  there  has  been 
overflow,  underflow,  or  neither.  Though  such  systems  are  rare,  they  have 
been  implemented  at  Toronto  [3)  and  at  Stanford  [2],  and  these  systems 
S'ake  programming  such  an  algorithm  as  QUADRATIC  far  more  satisfactory. 

Some  systems  merely  dismiss  a  user's  program  in  case  of  overflow. 

If  a  result  underflows,  many  systems  merely  set  the  result  to  0  and 
return  to  the  user's  program  without  any  indication.  Faced  with  systems 
like  these,  the  programmer  must  take  great  pains  to  insure  that  over-  or 
underflow  can  never  occur.  These  precautions  make  a  satisfactory  algorithm 
tedious  to  write,  lengthy  to  store,  and  slow  to  execute.  One  necessary 
subroutine  must  determine  the  exponents  of  a,  b,  c,  and  other  real 


numbers  local  to  the  procedure.  This  is  probably  most  gracefully  written 
in  machine  code. 

The  ability  to  achieve  the  prescribed  accuracy  (7)  of  0+1  units 
in  the  last  place  depends  on  the  detailed  properties  of  the  floating-point 
arithmetic  processor.  If  the  prescribed  accuracy  is  not  achievable^ 

condition  (7)  must  be  relaxed  as  necessary.  If  necessary,  one  could 

*  «■ 
also  make  the  set  F  smaller  by  changing  (3  ). 

As  a  simple  illustration  suppose  that  ?  =  10,  t  =  4,  m  «  -5**,  M  ■  U5  . 

Then  the  smallest  and  largest  positive  numbers  in  F  are 

1000  X  io‘^^  *  10"^^ 
and 

9999  X  10^^  1=  .9999  X  10^^  =  (1  -  10"^)  X  10^^  . 


Here  are  some  equations  that  may  give  trouble  for  this  system: 

(a)  10"^°z^  -  5  X  lo"^°z  +  6  X  10"^®  =  0; 

The  roots  are  2  and  3,  and  the  only  danger  is  that  undetected 
underflows  will  introduce  an  error. 

p  in 

(b)  z  +  10  z  -1  =  0: 

Use  of  the  standard  quadratic  formula  will  yield  0  for  the 
positive  root,  instead  of  the  correctly  rounded  value 


1.000  X  10"^°  . 

(c)  10 -  10^°  +1=0: 

70  -iiO 

The  roots  are  near  10  and  10  ;  use  of  the  quadratic 


formula  can  easily  cause  overflow  or  underflow. 
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(d)  2.BG^  -  2.864  z  +  0.?160  =  0: 

Here  0.5000  is  a  double  root.  Evaluating  the  quadratic 
formula 

2.864  +*\/(2.864)^  -  (4  X  2.864)  x  0.7160 
2  X  2.864 

in  single-precision  rounded  arithmetic  yields  0.5000  +  0.05477  i> 
with  an  error  of  over  547  units  in  the  last  placr  of  0.5000  . 

It  is  not  purely  academic  to  make  strict  demands  of  the  procedure 
QUADRATIC  .  Quadratic  equations  arising  in  the  course  of  solving  deter- 
minantal  equations  by  Muller's  method  L'6]  or  Daguerre's  method  [5],  parti¬ 
cularly  in  connection  with  large  matrices,  sometimes  have  one  of  the  roots 

* 

out  of  the  range  of  F  ,  and  yet  make  essential  use  of  the  root  in  the 

*  ,  . 

range  of  F  .  The  accuracy  requirement  (7)  is  perhaps  overstrict  for 
equations  with  nearly  double  roots. 

I  believe  the  specifications  to  be  very  reasonable  for  a  basic 
process  like  solving  a  quadratic  equation.  Nevertheless,  I  venture  to 
guess  that  not  more  than  five  quadratic  solvers  exist  anywhere  that  meet 
the  general  level  of  the  specifications.  Kahan  [4]  has  prepared  an 
algorithm  (in  Fortran  tv  for  the  7094-11  under  the  Toronto  version  of 
the  IBSYS  operating  system)  which  achieves  the  specifications  for  0=2, 
t  =  27,  m  =  -155,  M  =  100  .  The  error  never  exceeds  9A  3  +  lA) 
units  in  the  least-significant  place  of  the  root.  The  only  multiple- 
precision  operations  occur  in  the  computation  of  A  =  b^^  - 
followed  by  storage  of  a  single-precision  value  of  A  .  (Here  a^^,  b^^,  Cj^ 
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are  proportionai  to  a,  b,  c  .)  The  excess  time  for  the  double-precision 
computation  is  negligible  in  comparison  with  the  time  required  to  deal 
with  over-  and  underflow. 


It  is  obviously  relevant  to  ask  to  what  extent  the  various  computer 
algorithms  presented  at  this  SiTnposium  for  general  polynomials  make  pro¬ 
vision  for  over-  and  underflow,  and  also  what  accuracy  they  achieve. 

It  is  noteworthy  that  the  programming  of  QUADRATIC  depends  crucially 
on  the  arithmetic  properties  of  the  computing  system,  especially  on  its 
behavior  with  over-  and  underflow.  The  practical  numerical  analyst  with 
high  standards  is  thus  inextricably  involved  with  the  arithmetic  behavior 
of  his  digital  computer  hardware  and  accompanying  operating  systems. 
Unfortunately,  few  numerical  analysts  have  formulated  their  systems 
requirements  explicitly,  not  to  mention  communicating  them  effectively 
to  the  persons  who  design  hardware  and  software  systems.  With  existing 
computing  systems  a  numerical  analyst  faces  a  most  disagreeable  dilemma-- 
either  he  writes  less  than  satisfactory  algorithms,  or  he  undertakes  the 
never-ending  chore  of  writing  basic  software  systems  (and  perhaps  even 
rebuilds  the  arithmetic  unit) .  Professor  Kahan  has  generally  taken  the 
second  alternative. 


t 


References 

1.  George  E.  Forsythe  and  Cleve  B.  Moler,  COMPUTER  SOLUTION  OF 
LINEAR  ALGEBRAIC  SYSTEMS,  Prentice-Hall,  Englewood  Cliffs,  N.  J.,  1967} 
see  Sec.  20. 

2.  Michael  D.  Green,  "Coping  with  over/underflow  on  the  B5500, " 
multilithed  typescript.  Computer  Science  Dept.,  Stanford  Unl’.’crsity, 
Stanford  Calif.  9U305,  USA,  15  July  I966,  c.  50  pp. 

3.  W.  Kaha'i,  "709*+-II  system  support  for  numerical  analysis,"  draft. 
Department  of  Computer  Science,  University  of  Toronto,  Toronto,  Canada, 
Aug.  1966. 

1+.  W.  Kahan,  The  FORTRAN  IV  subroutine  QPRTC,  computer  library  of 
McLennan  Laboratories,  University  of  Toronto,  Toronto,  Canada,  c.  I966. 

5.  H.  J.  Maehly,  "Zur  iterativen  Auflosung  algebraischer 
Gleichungen, "  Z.  Angew.  Math.  Phys.,  ^  26C-263  (195^) • 

6.  David  E.  Muller,  "A  method  for  solving  algebraic  equations 
using  an  automatic  computer."  Math.  Tables  Other  Aids  Comput.,  10»  208- 
215  (1956). 


Computer  Science  Department 
Stanford  University 
Stanford,  California  9*^305 


Unclassified 


OONTtOLOATA-MD 


•IlMt. 


I.  eai«MATM«  ACTIVITV 

Computer  Science  Department 
Stanford  University 
Stanford,  California  9*^305 


a-  RSeONT  TITLt 


WHAT  IS  A  SATISFACTORY  QUADRATIC  BQUATICW  SOLVERT 


«.  DCaCM^IVI 

Manuscript  for  Publication  (Technical  Report) 


a.  AUTNOiirti  (tM*  MMM. 


TMCUMTV  CkMMeiea 


Unclasalfled 


George  E.  Forsythe 


a-  MPoar  dati 

Axigust  7,  1967 


Sa.  eONTHACT  OR  RRANT  NR. 

Nonr-225(37) 

k  RRRJBCT  *M. 


?•.  VeTAkMR.  RR  RARU  I  Va  MR. 

nine  I 


Sr  RMRMIATRR^  RSRRRT  IIUMRSm 


C8  7^ 


NR-Ol^U-211 


ta  AVAikAsii.iTv/i.iiirrATio«  MeTtesa 


RffRSRT  Nsm  (A 


none 


Releasable  without  limitations  on  dissemination. 


If  -  aURSkSMSMTARV  NOTSS 


la-  ASaTRACT 


imswaosiws  mmtarv  astvitv 

Office  of  Haval  Research 
Washington,  D.  C.  20360 


Cods  k32 


Following  the  suggestions  of  Professor  W.  Kahan,  the  author 
details  precise  requirements  for  a  satisfactory  computer  program  to 
solve  a  quadratic  equation  with  floating-point  coefficients.  The 
principal  practical  roblem  is  coping  with  overflow  and  underflow.  . 


Unclassified 


LINK  A 


ML!  I  «T 


quadratic  equation 

overflow 

iinderflow 

floating-point  computation 


t  omoiNA'tun  Acnvrrvt  sieat  uw  sm*  ms  as^ma 

aI  tlw  AAAlrAMAr.  AAtllAWAlHf,  ■TAAtAA^  PiflMAt  aID» 
Iaaaa  AAttfUf  m  9tkm  rngmlmUm  (oa^aeaf  AAtkAd  lAAAtag 

IlM  ftpAlt 

2a.  nratTIBCtMTYCLAMnCATIOMt  BaIwUmavac 
■U  AACAfitf  aIaaaMaaUaa  at  tlw  iAMil>  IaAIiaIa  wlmliw 
••Maatnatat  Oala"  to  UoIaSaA  MaMiv  to  M  b#  to  abaw* 
AAAA  aSIA  toswptton  AACHrMy  rAgtoAttAAA. 

aa.  OMXIPi  AAtoAiwto  SAATi^Aili^  to  yAcWAS  to  DaP  P» 
McUvA  SaOO,  10  aa4  AhaaA  Baacaa  lA<AAinAl  MaaaaI  Bator 
tha  WAMA  aaabar.  AIaa.  vlaia  aoeltoabto.  aIwa  Uwt  tpttoAal 
awrirlAja  kava  baaa  aaaO  far  Qr^  S  ami  Otaao  4  aa  atUwr* 


ritr  daaalflcattoA.  Aatog  ati 


a.  RBPOBT  TITLBi  Baiar  tha  AAAplata  rapait  Utto  to  all 
aapltol  tottara.  Tltlaa  to  aU  aaaaa  AlwAie  ha  AAalAAAtflaAi 

U  A  AWAAil^l  titto  CAWWt  bA  AAlACtoO  VltlWAt  CltoAtWC^ 
tto%  Ahaar  Ulla  dAAAlWcAttoA  to  all  caqIIaIa  to  yawatbAAto 
baawAtotaly  Mlaarias  tha  tlUa. 

4.  PBKRlPnVB  NOTBP  U  aoprapriata.  aatar  tha  typa  aC 
•apart,  a.g.,  iwarha,  prapraaa,  AMwaara,  aawial,  ar  ftoat 
Qlaa  tha  taduAlaa  Oataa  artwa  a  apaalfic  rapaitlag  parto4  to 


S.  AVrHOR(S)i  Batar  tha  aaaH(a)  at  aathaf(a)  aa  ahaara  aa 
ar  to  tha  rapait  Batar  laat  iwaa,  drat  aaaw,  wMBa  tottlal. 

If  aUlMary,  ahaar  raiOi  amd  hraiwh  af  aaralcAb  Tha  aajw  at 
tha  prtoel^  authar  to  aa  abaatota  ■lahAAAi  ratBlraaaat 

A.  RBPORT  OATD  Batar  tha  4ata  at  tha  rap  art  aa  4ay, 
aiaacK  faar,  ar  awath,  jraw^  If  'baa  aaa  4ala  appaara 

oa  tha  rapait  aaa  Saw  af  pahltcattoib 

7a.  TOTAL  NUIBBII  OP  PAOBSt  Tha  taCd  p^  caaal 
ahauM  fatlew  aanaal  paptoatlaa  pracaOwraa,  tab,  aatar  tha 
iMatoar  af  pagaa  caatalalag  toforawtta» 

7A.  NIUBER  OP  RBPBRBNCBP  Batar  tha  fatal  aaabar  at 
rafaraiwaa  cita4  la  tha  rapait, 

•a  CONTRACT  OR  GRANT  NUM8BR:  If  apprapriata,  aatar 
tha  appllcabla  auatoar  of  tha  eoatract  ar  ^aal  uitoar  toilch 
tha  rapen  waa  wlttaA 

•A.  to,  to  $d.  PROJECT  NUMSER:  Batar  tha  apprapriata 
aUlltary  dapartowat  IdaatlflcaUaa,  aaek  aa  prajact  aoarbar, 
auhproiact  auabar,  ayataw  auatocra,  taak  aaatoar,  ate. 

Oa.  ORIGINATOR’S  REPORT  NUIIBER(S):  Batar  tha  aftt* 
clal  raporl  ouabar  by  which  tha  docuataat  will  ha  toaallflad 
and  controlled  by  tha  orlalaattoii  actletty.  Thia  aaabar  awat 
ha  ualqua  to  Ihla  raport 

9b.  OTHER  REPORT  NUMBEfKS):  II  tha  raport  haa  baan 
aaalpnad  any  oihrt  raport  nuaibara  f  allbar  by  Iba  ertglnamt 
ar  by  Ibo  aponaor),  alao  aniar  thla  aiaabarfaX 

10.  AVAILADILITY/LIUITATION  NOTICES;  Batar  aay  lira- 
Itotlona  on  lurihci  dlaaaatoMttoa  af  tha  raport,  othar  thaa  Umm 


DD  .'..r..  1473  (BACK) 


(I)  "Qaallflad  laeaaatarA  awy  abtato  aaptoa  ad  thto 
rapaat  Raai  DOC" 


rapartiy  DDC  to  aat  aatharlaasl""*"*****  ** 

(JD  "U  I.  Oaaaraaaat  ataadaa  way  ahtato  aaptoa  af 
thla  rapart  dhaetly  Aaai  DOC  Othar  faaimad  DOC 


(4)  "O.  to  adUtaryapaacitotoay  ahtato  aaptoa  af  thla 
lapart  Rlractly  ha«  DDC  Odhar  paalUlaS  aaara 


(S)  "All  dlatrtoatlaa  af  thla  rapart  la  c.  .aO  Qsd> 
WaS  DOC  aaara  ahall  raaaaat  thraiwh 


If  tha  rapait  haa  baaa  faralahad  to  tha  OfRea  af  Taahalaal 
•aratoaa,  Dapartaiaat  af  Caaaaaraa,  far  aala  to  tha  paMte,  ladb 
aata  thla  faat  aad  ardar  tha  prlaab  If  kaaww 

IL  fUPPLBMBNTARY  NOTBR  Uaa  far  addHIaaal  a^lMto 


IX  SPONKMUNO  MILITARY  ACnvITT)  Batar  tha  aaM  af 
tha  dtp  Alfa  atal  prajact  afSea  ar  labaratary  apaapartaa  Cmp 
Irtg  lat)  tha  raaaarch  aad  daaalayaat  laciada  adtoaaai 


Irtg  tai)  tha  raaaarch  aad  daralayaat  laciato  adtoaaai 

IS.  ABSTRACT  I  Batar  aa  abatiaet  ylviap  a  brtof  aad  laetaal 
aaaaaary  of  iha  dacaaaat  ladtoatlra  of  tha  rapart  atraa  thaagh 
It  tiay  atoa  appaar  olaaadwia  la  Mw  ba^  of  tha  Mahalaol  ra* 
part  K  addlUaiwI  apaca  la  roeoirad,  a  aaatlaaatlaa  ahoal  aha 
ha  attachad. 


It  to  hlAly  daalrabla  that  tha  ahatraet  af  tlaaoUlad  wpaita 
ha  aaetoaoinad.  Bach  parapiaph  of  dw  abotioet  ahall  aad  artih 
aa  todtoatlaa  af  tha  aUUlaiy  atHaarlty  etoaalfloattoa  af  tha  la- 
fanoaUoa  to  tha  yiayyh,  lapraaaalad  aa  frp>.  fS>.  fCJ,  or  fWj. 


Thara  la  aa  liwilaUaa  aa  dm  laagth  af  tha 
atror,  tha  aopgaatad  loagib  la  ham  ISO  ta  22$  a 


ahatraet  Maw- 


14.  KBY  WORDS:  Koy  ararda  arc  tochatoally  awaatapM  taiM 
ar  abort  pbiaaaa  that  cbaractartoa  a  rapart  aad  atap  ha  aaad  aa 
todaa  aatriaa  far  catalegtaf  the  raport.  Bay  worda  aatat  ha 
aalactod  ao  dwt  aa  aacoriiy  cloaolOcatlM  la  raaalrad  Idaatl- 
ftoro.  Back  aa  aqalpaMol  aMdol  doalipatloa,  tiada  aawa,  toUltary 
projMt  coda  Mwo.  gaotiaphlc  location,  aay  bo  aaad  aa  hoy 
wordo  hat  will  bo  followod  by  an  todIcaUoa  of  tochalcal  caa- 
toat.  Tbo  aaalpaawai  af  Uaka,  lalaa,  aad  wol^a  to  aptlaaal. 


Unclassi  Tif-d 

Seciirity  ClasaiOcatic 


